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We report the first application of critical cluster techniques to the Mott metal-insulator transition 
in vanadium dioxide. We show that the geometric properties of the metallic and insulating puddles 
observed by scanning near-field infrared microscopy are consistent with the system passing near 
criticality of the random field Ising model as temperature is varied. The resulting large barriers to 
equilibrium may be the source of the unusually robust hysteresis phenomena associated with the 
metal-insulator transition in this system. 


The Mott metal-insulator transition in VO 2 has the po¬ 
tential to produce disruptive technologies, such as mem- 
ristors, memory capacitors, ultrafast switches, and possi¬ 
bly even neuromorphic circuits [1]. Yet despite decades of 
research into the Mott metal-insulator transition in VO 2 , 
the nature of the phase transition, and in particular the 
broad hysteretic behavior accompanying it, is not yet un¬ 
derstood. The Mott metal-insulator transition has some 
universal features of the liquid-gas transition [2, 3j: the 
transition occurs through a first order phase transition 
line broadened by conductivity hysteresis, which termi¬ 
nates at a classical critical point associated with general 
Ising universality 12 6 ]. In addition, in VO 2 a percolation 
model has been proposed through the study of resistance 
avalanches [7], while in V 2 O 3 the clean 3D Ising model 
has been implicated [2]. By applying quantitative crit¬ 
ical cluster techniques to study the criticality, we show 
that a key ingredient missing from prior treatments is 
the prominent role of disorder, with important implica¬ 
tions for the robust hysteresis effects associated with the 
metal-insulator transition. 

As shown in Fig. [lj scanning near-field infrared mi¬ 
croscopy (SNIM) on VO 2 0 reveals the complex pat¬ 
tern formation associated with the transition from the 
low temperature insulating phase to the high tempera¬ 
ture metallic phase. We apply recently developed clus¬ 
ter techniques [9] to the observed multiscale patterns of 
inhomogeneous local conductivity, revealing critical ex¬ 
ponents such as the fisher exponent r, volume fractal 
dimension d v , and hull fractal dimension dy, both for 
clusters and avalanches. This method has the advantage 
that it can access the true universality class of the phase 
transition, unlike methods based on macroscopic conduc¬ 
tivity in which the exponents are affected by the contrast 
ratio between metallic and insulating regions :3]. Our 
analysis is the first to quantitatively incorporate disorder 
into the study of the criticality, allowing us to determine 
the relative importance of disorder and interactions, and 
identify the dominant type of disorder. 



FIG. 1. Scanning Near-field Infrared Microscopy on VO 2 as 
temperature is increased through the Mott metal-insulator 
transition regime, over the same 4 fim x 4 /im sample area. 
These figures are Ising-mapped from the original SNIM im¬ 
ages in Ref. [8] with threshold scattering amplitude ath = 2.5. 
The metallic regions, colored green, give higher near-field 
scattering amplitude, a > a t h , compared with the insulating 
regions with a < ath , colored blue. 


Near a critical point, the correlation length grows to 
become the dominant length scale, and it is possible to 
map the real physical system to a coarse-grained model 
with the same universal features. In the original SNIM 
data images 18], we assign Ising variable a = 1 (metallic) 
or —1 (insulating) to each coarse-grained region @- 0 |. 
The threshold amplitude a t h for identifying the metallic 
regions (those with scattering amplitude a > a t h,) and 
insulating regions (with a < a t h) is about 2.5 la El, 
which we use throughout this Letter. We furthermore 
incorporate disorder into the model: 
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where the sum runs over the coarse-grained regions (sites) 
consisting of a cubic lattice, chosen with spacing at least 
as small as the resolution of the images to be studied. 
The tendency for neighboring regions to be of like char¬ 
acter is modeled as a nearest neighbor ferromagnetic in¬ 
teraction J > 0. Because the data considered is that of a 
thin film, the sum over Ising variables in the plane of the 
film (denoted by ||) extends to infinity, but the sum over 
Ising variables perpendicular to the film surface (denoted 
by _L) is finite, confined by the film thickness L z . De¬ 
pending on the size of the correlation volume compared 
to the size of the system, the sample may display two- 
dimensional (2D) or three-dimensional (3D) critical be¬ 
havior near criticality. 

At the order parameter level, there are two broad 
classes of disorder: local energy density disorder (which 
we incorporate as random bond disorder), and random 
field disorder & Random bond disorder is included 
through the term SJij, and hi represents random field 
disorder, which is chosen from a gaussian probability dis¬ 
tribution centered about zero, with variance R. R is often 
called the disorder parameter or just disorder. The field h 
represents a generalized external field which couples with 
the local Ising variables. Figure [ 2 ] shows the schematic 
phase diagrams associated with this general Ising model 
of Eqn. [U emcompassing several universality classes both 
in 2D and 3D. Note that random field disorder is always 
relevant in the renormalization group sense, and in fact if 
both random bond and random field disorder are present, 
the associated stable fixed points are those of the random 
field disorder. 

Using our mapping to Ising variables (Fig. [lj), we track 
the geometric clusters, defined as connected sets of near¬ 
est neighbor sites (pixels) with the same color. We then 
use the statistics of the sizes and shapes of these geomet¬ 
ric clusters to identify the cause of the complex pattern 
formation. In comparing the spatial complexity revealed 
in the SNIM data to theory, there are 7 fixed points to 
consider, as shown in Fig. [2] In the two-dimensional case, 
these consist of the clean Ising model (C- 2 D, i.e. in the 
absence of any material disorder), uncorrelated percola¬ 
tion (P-2D), and the random field Ising model (RF-2D). 
Note that random bond disorder is irrelevant in the renor¬ 
malization group sense in 2D, and the phase transition is 
therefore governed by clean Ising model exponents set by 
the C- 2 D fixed point. Although the RF- 2 D fixed point is 
unstable, its critical behavior can be observed for weak 
enough finite disorder. In the three-dimensional case, 
the possible fixed points are the clean Ising model (C- 
3D), uncorrelated percolation (P-3D), the random bond 
Ising model (RB-3D), and the random field Ising model 
(RF-3D). For the exponents extraction, in order to re¬ 
duce noise due to the finite field of view (FOV), we apply 
the logarithmic binning method throughout this Letter, 
which is a standard technique for analyzing power law be¬ 
havior 0 . The value of critical exponents are obtained 




FIG. 2. Schematic equilibrium phase diagrams and fixed 
points of Eqn. [I] Two classes of disorder: (a) Random field 
disorder (with or without random bond disorder); (b) Ran¬ 
dom bond disorder in the absence of random field disorder. 
Solid regions denote the ordered phase, from 2D (yellow re¬ 
gion, applicable when the correlation length £ L z ) to 3D 
(green region, applicable for £ large but less than L z ). Solid 
green lines represent continuous phase transitions. The green 
arrow (flow) on each phase transition line points to the solid 
green circle representing the fixed point controlling the long 
distance (universal) behavior of that phase transition. 


by taking the discrete logarithmic derivative (DLD) 0. 

From the SNIM images, we extract three critical ex¬ 
ponents, which are the Fisher exponent r, the volume 
fractal dimension d v , and the hull fractal dimension dh, 
as entailed by the self-similarity of the geometric clusters 
near certain critical points, r characterizes the cluster- 
size distribution D(s), which is the histogram of cluster 
sizes s, and scales as D(s) oc s~ T . d v and dh characterize 
the fractal nature of cluster sizes s and hulls h. For the 
cluster sizes we have s oc R^ v , where R s is the radius of 
gyration of the cluster [14]. For the cluster surfaces, we 
have h oc R^ h . Rh here refers to the radius of gyration 
of all the sites enclosed by the hull, including any sub¬ 
clusters inside. Throughout this letter, we analyze the 
power law behaviors using only the internal clusters, i.e., 
clusters which do not intersect the boundary of the FOV, 
in order to mitigate finite FOV effect to the extraction of 
exponents. Although from our analysis d v and dh do not 
appear to be affected by this effect, estimates of r in a 
finite FOV are always skewed to lower values because of 
a pronounced bump irqthe scaling function of the clus¬ 
ter size distribution [10, [H, 16]. Within a cutoff in the 
decades of scaling, D(s) of internal clusters would share 
the same power law as the full system [l 6 |, therefore r 
could be extracted more accurately. 

Figure [3] shows the extraction of r, d v , and dh by ap¬ 
plying our cluster techniques to the finite size SNIM im¬ 
age data through the Mott transition using internal clus¬ 
ters. The main panels show explicitly the power law fits 
of these critical cluster exponents at three intermediate 
temperatures (which have enough clusters for good statis¬ 
tics) with threshold scattering amplitude a t h = 2.5. For 
r, two decades of scaling are evident, and the remaining 
points fall off the scaling regime due to the use of internal 
clusters 0 ; and for the fractal dimensions, with the first 
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FIG. 3. Power law fits of internal clusters from the VO 2 SNIM 
imaging datasets with threshold scattering amplitude ath — 
2.5 for (a) D(s) oc s _T , (b) R s oc s 1 /^, and (c) Rh oc h 1 / dh , at 
three intermediate temperatures T — 342. 6K, T = 342.8 K, 
and T m. 343 K. The insets show r, l/d v , and 1/dh (with 
error bars) extracted using the same method as a function of 
ath at T=342.8K. 


bin (which contains s = 1 clusters) excluded 10]. robust 
power law scaling extends over multiple decades, encom¬ 
passing the entire FOV. The robust power law behaviors 
as well as the fact that the values of the cluster exponents 
are the same within error bars for different intermediate 
temperatures corroborates the idea that the system is 
near some critical end point of the Mott transition. The 
insets of Fig.[3]show our extracted r, d V: and dh using dif¬ 
ferent ath at the representative temperature T = 342. 8K 
which is closest to criticality, and within error bars these 
critical exponents are robust against changes in threshold 
amplitude within a broad stable region around a t h = 2.5, 
consistent with that given by the definition in the exper¬ 
iment paper BEcf This independence of results with 
respect to microscopic details also reflects universal be¬ 
havior near criticality. Using a t h = 2.5 and T=342.8K 
as representative conditions, we have r = 1.92 ±0.11, 
d v = 1.93 =b 0.07, and dh = 1.21 zb 0.03 for the spatially 
complex clusters near the critical end point of the Mott 
transition in VO 2 . 

Before comparing the numerical values of the critical 
exponents revealed by our analysis of the data to the 
known theoretical values of each fixed point, the C-3D 
and RB-3D fixed points can be ruled out a priori. Since 
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FIG. 4. Exponent comparison charts, (a) Critical exponent 
r. (b) Critical exponents d v (blue) and dh (green). The hor¬ 
izontal fines are our extracted exponents from the SNIM im¬ 
ages, with the shaded regions being their error bars. The 
circles represent theoretical values HE} for the fixed points 
of Eqn. [lj When comparing with 3D models, we have to use 
the effective values of d* and d ^ corresponding to taking 2D 
cross sections of the clusters embedded in 3D, i.e. d* = 2d v /3 
and d^ = dh /2 S M- These effective values are represented 
by the open circles. 


geometric clusters are not fractal at these two fixed points 
17-19], they cannot be the cause of the robust power 


law behavior observed in Figs.[3fb) and[3jc). Five fixed 
points remain to be considered (C-2D, RF-2D, RF-3D, 
P-2D, and P-3D), and the geometric clusters do exhibit 
a fractal nature at these fixed points Bill ELM. In 
the slab geometry considered here (and also for layered 
materials), it is also necessary to consider possible di¬ 
mensional crossovers. Within this set of candidate fixed 
points, there are two possible dimensional crossovers. In 
the presence of random field disorder, the dimensional 
crossover is between RF-2D and RF-3D. (For a slab ge¬ 
ometry, exponents will drift from the 3D fixed point to¬ 
ward the associated 2D fixed point, once the correlation 
length in the c-axis direction begins to exceed the slab 
thickness.) The other possible dimensional crossover is 
between C-2D and the 3D percolation points T p < T c in 
the clean and random bond Ising models. Although the 
exponents we require have not been explicitly reported 
in the literature, the expectation in the literature is that 
the exponents should follow those of uncorrelated perco- 
lation, P-3D [19,El,E2]. 


As shown in Fig. 01 the values of r and d v resulting 
from our analysis of the data are fairly close to all fixed 
points other than P-3D. Therefore, uncorrelated 3D per¬ 
colation can be ruled out as an origin of the power law 
behavior of the statistics and geometry of the metal and 
insulator islands. Note, however, that the theoretical val¬ 
ues of the remaining four fixed points are all fairly close in 
value for these two exponents, making it hard to distin¬ 
guish among the remaining fixed points using only these 
two measures. On the other hand, dh varies significantly 
for the various fixed points. We see immediately from the 
comparison chart that the theoretical values of dh at P- 
2D and RF-2D are inconsistent with the data. Whereas 
the closest match for this measure is P-3D, that is incon- 
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FIG. 5. Power law fits of internal avalanches from the VO 2 
SNIM imaging datasets with threshold scattering amplitude 
a t h = 2.5 for (a) D(s) oc s~ T , (b) R s oc s 1 ^ d ' v , and (c) 
Rh oc h}^ dh , at three intermediate intervals with AT = 0.2 K 
through the Mott transition 


sistent with the other two critical exponents, and must 
remain ruled out. 

Only two candidate fixed points remain (C-2D and RF- 
3D), but in Fig. 0(b), both show about a 13% discrepancy 
between the data and the theoretical model for the value 
of dh , which is significantly larger than the typical error of 
about 3% of the exponents extracted from the data. We 
turn then to the possibility of dimensional crossovers. In 
the case of the clean or random bond Ising model, power 
law behavior of geometric clusters will drift from P-3D to 
C-2D as larger length scales are observed. However, the 
value of d v is much closer to C-2D than to P-3D, consis¬ 
tent with the geometric clusters being near their 2D limit, 
whereas the value of dh is much closer to P-3D than to 
C-2D, which would indicate that the geometric clusters 
are near their 3D limit. This inconsistency strongly ar¬ 
gues against the P-3D to C-2D dimensional crossover as 
the origin of the power law behavior. 

The other candidate dimensional crossover is from RF- 
3D to RF-2D. The exponent comparisons show no incon¬ 
sistency with this explanation, since within error bars 
r and d v match the whole RF-2D to RF-3D crossover 
regime, and dh is also consistent with this dimensional 
crossover. Therefore, the random field Ising universal¬ 
ity class best describes the critical behavior of the Mott 
transition in the VO 2 thin film, revealing the key role 
played by disorder effects in these systems. Observations 
at longer length scales can be used to test this hypothe¬ 
sis: if our identification is correct, then the values of all 
exponents should drift toward RF-2D with larger fields 
of view in a film geometry such as the present one. 

In order to further test whether the spatial complex¬ 


ity is driven by random field effects, we also analyze the 
avalanches which, in the context of the VO 2 experiment, 
are defined as the difference between Ising maps between 
two neighboring temperatures. We study the critical 
avalanche exponents arising from the self-similarity of 
avalanches and are directly extractable from the exper¬ 
imental image datasets. Figure [5] shows the power law 
fits of the internal avalanches with threshold a t h = 2.5 
for the data-extracted r*, d*, and dj£, using the same 
techniques as for the clusters in Fig. [3j The star symbol 
is used to denote that these exponents are for avalanches. 
Similar to their cluster counterparts, the cut-off for in¬ 
ternal avalanche scaling for r* is about 1.5 to 2 decades 
of scaling, and d* and d^ show multiple decades of power 
law scaling behavior. For different intervals, the ex¬ 
tracted critical avalanche exponents are consistent with 
each other within error bars. This data-extracted robust 
universal power law scaling of avalanches provides further 
evidence that the random field universality is underlying 
the complex pattern formation at the Mott transition in 
the VO 2 thin film, since non-trivial scaling behavior of 
avalanches is characteristic of random field physics. 

Since the avalanches are taken from the finite interval 
between scanned images, the data-extracted t* in our 
context characterizes the scaling of the field integrated 
avalanche size distribution Di nt (s ) near criticality. The 
corresponding exponent is reported to be 2.03 =b 0.03 for 
RF-3D [HHlQ and 1.3 ±0.1 @ for RF-2D, and the 
corresponding data-extracted one from Fig. [5(a) lies in 
between, consistent with the conclusion of a dimensional 
crossover from RF-3D to RF-2D. With d* and d £ only re¬ 
ported for RF-3D [24], our data-extracted fractal dimen¬ 
sions Fig. [5(b) and (c) could not be directly compared 
with the theoretical ones. However, by comparison with 
Fig. [3(b) and (c), within error bars d v and dh extracted 
from avalanches have the same values as those from clus¬ 
ters, consistent with Ref. 0 , which show that clusters 
and avalanches have the same exponents in random field 
Ising model (RFIM). 

The dominance of random field disorder is quite rea¬ 
sonable in a real physical system, arising from defects 
and impurities in the crystal structure; randomness in 
grain size, orientation, and boundary structure; and im¬ 
perfections in the substrate. In particular, it has been 
shown that the temperature at which grains transition 
depends on grain size [26]. In addition, we have previ¬ 
ously shown that metallicity nucleates first near the grain 
boundaries (see Fig. 6 of Ref [27]). Since physical tem¬ 
perature maps to effective field in our model, these effects 
map to random field disorder in the model. (See SI for a 
more thorough discussion.) One prediction of this line of 
reasoning is that grain size, since it affects the random 
field strength in the model, is expected to correlate with 
the hysteresis width. 

Several other characteristics of the data serve to cor¬ 
roborate the random field hypothesis. For example, as 
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the Mott transition proceeds, metallic puddle formation 
occurs via nucleation and not exclusively by front prop¬ 
agation [28|, as can be seen in Fig. [T] In addition, as 
temperature increases through the transition (analogous 
here to sweeping the generalized field h of Eqn. [I]), do¬ 
mains change from a = — 1 to cr = +1, but never re¬ 
vert. This no-passing rule is consistent with the random 
field universality class 13, 0 0, indicating that the 
effects of quenched disorder dominate over thermal fluc¬ 
tuations. The dominance of quenched disorder over ther¬ 
mal effects is further corroborated by the reproducibility 
of the SNIM images, in that repeated near-field scans 
in the insulator-to-metal transition regime over the same 
sample area and at a fixed temperature show nearly iden- 
tical patterns of metallic puddles in the insulating host 
@, 0 . 

Finally, perhaps the most tell-tale sign of random field 
behavior in the data is the large width of the hystere¬ 
sis in temperature, about 7.5 K for temperature sweeps 


taking several minutes in a thin film |3jJ, [32[. The ex¬ 
treme critical slowing down of the random field case 
means that barriers to equilibration grow much more 
severely as criticality is approached than the usual “crit¬ 
ical slowing down” would predict. Hysteresis typically 
becomes long-lived near the critical endpoint of a first 
order phase transition, with the relaxation time t re i di¬ 
verging as a power law near the critical endpoint (33| . 
t re i oc oc \g — g c \~ uz : where g represents the relevant 
variable, whether temperature T or disorder strength 
R. £ is the correlation length, v is the exponent of 
the correlation length, and 2 is the dynamical exponent. 
Rather than a mere power law, barriers to equilibration 
grow exponentially near the RFIM critical endpoint [34|, 
ex PK 0 ] oc exp[l/|# ~ 9c\ u0 ], where 0 is the ”vio¬ 
lation of hyperscaling” exponent. These exponentially 
large barriers to equilibration are consistent with the 
large width of hysteresis loops evident in this Mott metal- 
insulator transition. In addition, the lar ge c ritical region 
typically associated with RFIM physics [15| is consistent 
with the wide range of parameters over which the broad 
hysteresis is observed. 

In conclusion, by applying newly developed cluster 
techniques to SNIM data on VO 2 , we have shown us¬ 
ing three different critical exponents that the critical end 
point of the Mott transition is in the universality class 
of the random field Ising model. This finding reveals 
a delicate interplay between interaction and disorder in 
these systems. The random field Ising universality class 
has also recently been shown to account for the spatial 
structure of the locally oriented domains observed via 
scanning tunneling microscopy on cuprate superconduc¬ 
tors (3, 1 


This further emphasizes the important role 
played by disorder in strongly correlated electron sys¬ 
tems, and indicates that there may be universality to 
the spatial complexity [36, observed in a wide vari¬ 
ety of strongly correlated electron systems. The cluster 


techniques employed here can readily be applied to 2D 
images in the context of other materials and microscopy 
techniques for the study of critical behavior. 
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